Magnetic order and paramagnetic phases in the quantum J1-J2-J3 honeycomb model 
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Recent work shows that a quantum spin liquid can arise in reahstic fermionic models on a honey- 
comb lattice. We study the quantum spin-1/2 Heisenberg honeycomb model, considering couplings 
Ji, J2, and J 3 up to third nearest neighbors. We use an unbiased pseudofermion functional renor- 
malization group method to compute the magnetic susceptibility and determine the ordered and 
disordered states of the model. Aside from antiferromagnetic, coUinear, and spiral order domains, 
we find a large paramagnetic region at intermediate J2 coupling. For larger J2 within this domain, 
we find a strong tendency to staggered dimer ordering, while the remaining paramagnetic regime 
for low J2 shows only weak plaquet and staggered dimer response. We suggest this regime to be a 
promising region to look for quantum spin liquid states when charge fluctuations would be included. 

PACS numbers: 75.10.Kt, 75.10.Jm 



The search for a quantum spin liquid phase in nature 
ever since has been a complicated task not only from ex- 
periment, but also from theory pi, This stems from 
the fact that the properties of a quantum spin liquid 
are very peculiar: quantum fluctuations have to form a 
many-body singlet state without long-range correlations 
of any kind of operator. Since the notion of frustration in 
quantum spin systems which is the main resource of fluc- 
tuations to accomplish a magnetically disordered quan- 
tum state at zero temperature, studies of a plethora of 
spin Hamiltonians on different lattices tell us that even 
most frustrated quantum systems tend to establish some 
sort of long-range correlations such as seen in a valence 
bond solid phase: while local spin operator correlations 
rapidly fall off there, long-range dimer-dimer correlations 
spoil the phenomenology of a quantum spin liquid. 

The remarkable numerical studies by Meng et al. 
are a fortunate exception. In their work, they report 
on the first unambiguous discovery of a genuine spin liq- 
uid phase from a generic microscopic model. They con- 
sider the Hubbard model on the honeycomb lattice by 
Monte Carlo methods and find a spin-gapped phase at 
U/t = 4.3 which shows no long-range correlations of any 
kind, neither charge density wave, superconductivity or 
even spin solid-type correlations such as that of a valence- 
bond crystal formation. Moreover, the study finds a clear 
excitation gap and no symmetry breaking of any lattice 
symmetry or parity P and time-reversal T, which already 
excludes chiral spin liquids and algebraic spin liquids. 

One path of further understanding the magnetic quan- 
tum phases on the honeycomb lattice is the develop- 
ment of effective descriptions for the Hubbard model 
itself such as gauge theory i4j and slave boson mean 
field theory methods O E]. Another direction, how- 
ever, which we pursue in this Letter is to analyze the 
Gutzwiller-projected Hubbard model on the honeycomb 



lattice. Specifically, we consider Heisenberg spin cou- 
plings up to third nearest neighbors labeled as the Ji- 
J2-J3 model. The motivation for this is two-fold. First, 
the J1-J2-J3 model projects onto the square root fraction 
of the Hilbert space of single site occupancy where only 
spin modes are present, which enables us to analyze the 
model through methods designed for this setup. Second, 
in reverse, a detailed understanding of the J1-J2-J3 phase 
diagram will eventually help to identify which aspects of 
possible quantum phases may be explained through spin 
fluctuations only and which may necessitate the effect of 
charge fluctuations. 

We consider the Hamiltonian 

^HCM = -^1 ^ ^i^j + -^2^ SiSj ^ Js^^^ SiSj, (1) 

ihj) {{hj)) {{{hj))) 

where the sums extend over nearest, second nearest and 
third nearest neighbors, respectively (see Fig. [T]). This 
model is expected to describe a class of magnetic ma- 
terials with a honeycomb lattice, one example being 
P — CU2V2O7 [7 . J2 and J3 are given in units of Ji. 
The solution of the classical J1-J2-J3 model has been 
known for a long time [8j [9]. For small J2 the system 
is Neel-ordered, which is commensurate with the honey- 
comb lattice and preserves the sublattice 120° degrees 
rotational symmetry. For sufficiently low J3 and beyond 
a threshold of J2, the system resides in a spiral phase. 
For high J2 and J3, it is energetically preferred to order 
in a collinear phase where spins along zigzag chains align 
ferromagnetically while neighboring zigzag chains exhibit 
antiparallel spin orientation (there are three degenerate 
collinear configurations). 

In order to obtain an adequate quantum phase diagram 
of the model defined in ([T]), there are not many suitable 
methods available, some of which predominantly focused 
on the J1-J2 line. Exact diagonalization (ED) studies 
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FIG. 1. (Color online) Phase diagram of the J1-J2-J3 honey- 
comb model via PFFRG (with exchange couplings depicted 
in the inset). The points represent parameter settings which 
we have computed. In the depicted J2-J3 range we find AFM 
Neel order (red circles), collinear order (C-AFM, green trian- 
gles), spiral order (blue squares) and a paramagnetic phase 
(open circles). The spiral order phase partly shows incom- 
mensurability shifts (open blue squares) from the spiral phase 
for dominant J2 (see also Fig. [2]). 



are very helpful, as arbitrary m-point correlation func- 
tions can be computed in principle. However, except of 
the valence bond crystal domain where the dimer Hilbert 
space projection [10 is justified [11 , ED cannot reach 
sufficient system sizes to adequately determine all phase 
regimes pTHl3] . Linear spin wave theory gives a qualita- 
tive tendency where the quantum corrections lead to, but 
is still strongly biased towards the classical limit p!2UT4] . 
Similar shortcomings apply to Schwinger boson mean 
field theory [13, 15j. A promising direction has recently 
been given by variational Monte Carlo (VMC) schemes, 
where energies of the antiferromagnetic (AFM) ordered 
trial state as well as several spin liquid (solid) candidates 
have been compared for J3 = 0, < J2 < 0.5 [16]. There, 
a transition point of J2 = 0.08 is found from the AFM 
ordered phase to the spin liquid phase which breaks no 
rotational invariance. It is followed by a transition to a 
dimerized spin solid phase at J2 = 0.3 which breaks rota- 
tional invariance of the lattice. It is hence likely that the 
disordered phase itself, already for the J1-J2 line and as 
such definitely for the full phase diagram with finite J3, 
contains different types of paramagnetic phases, which 
we will investigate in the following. 

We employ the pseudo-fermion functional renormaliza- 
tion group (PFFRG) [ITj il8j to identify magnetically or- 



dered and paramagnetic phases and compute the ground- 
state magnetic susceptibility to determine the magnetic 
order parameter in the ordered phases as well as the 
dominant fluctuation profile in the paramagnetic phases. 
Our starting point is the pseudofermion representation of 

spin-1/2 operators = 1/2 E«/3 /l^a/?//^' (^^^/^ =t4, 
/i = x^y^z) with the fermionic operators and and 
the Pauli-matrices cr^ . This representation enables us 
to apply Wick's theorem leading to standard Feynman 
many-body techniques. Quantum spin models are inher- 
ently strongly coupled models, requiring an infinite self- 
consistent resummation theory. In this context FRG [18]- 
|24 provides a systematic summation in different inter- 
action channels by generating equations for the evolu- 
tion of all one-particle irreducible m-particle vertex func- 
tions under the flow of an IR frequency cutoff A. In 
order to reduce the infinite hierarchy of equations to a 
closed set, we restrict the computation to the full set 
of one-loop parquet diagrams and their vertex correc- 
tions up to infinite order, and additionally include cer- 
tain two-loop contributions that are essential to induce a 
self-consistent resummation procedure exceeding the per- 
turbative limit [T8[ [25] . The parquet diagrams include 
graphs that favor magnetic order and those that favor 
disorder tendencies such that in total the method pro- 
vides an adequate treatment of order and disorder fluctu- 
ations [iSl [26l [27] . From the two-particle scattering ver- 
tex we obtain the spin-spin correlation function (or spin 
susceptibility), which is the central outcome of our PF- 
FRG. A magnetic ordering instability is initially signalled 
by a strong rise of the susceptibility associated with this 
order at some finite scale of A. The onset of spontaneous 
long-range order is signalled by a sudden breakdown of 
the smooth flow (as shown for J2 = 0.1 in Fig. |2|. In the 
momentum-resolved magnetic susceptibility x(k), where 
k is defined along kx^y components in the extended Bril- 
louin zone (BZ), the different magnetic ordering patterns 
manifest as peaks depicted in Fig. [2j Note that due to 
the two-atomic unit cell, the spin susceptibility does not 
have the periodicity of the first BZ but rather of the ex- 
tended (second) BZ. Adding small dimer-field perturba- 
tions to the Hamiltonian, we are able to calculate dimer 
responses for a given dimer configuration (Fig. [2|. For 
the present study our PFFRG algorithm internally deals 
with spin-spin correlations up to a length of 9 lattice 
constants corresponding to a correlation area of 181 sites 
which provides an adequate description of the model. 

Fig. [ijshows the quantum phase diagram of the J1-J2- 
J3 model ([1]). For dominant Ji the system displays AFM 
order which persists longer against J2 for finite J3 as J3 
cooperates with Ji. Increasing J2 (for not too large J3) 
we observe a melting of the order and the appearance of 
a rather large paramagnetic region. Above J2 ~ 0.6 the 
system is characterized by presumably weak magnetic or- 
der and very small ordering instability scales which are 
hard to resolve numerically. However, as we enter this re- 
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FIG. 2. (Color online) The Heisenberg honeycomb model along the J3 = axis. Upper line: J2-sweep of the static k-space 
resolved susceptibility. In magnetic phases (J3 = 0.1) the susceptibility is depicted at a cutoff scale just before the instability 
breakdown, otherwise the correlations are derived from the final flow values at A = 0. Susceptibilities are always given in 
units of j^. Lower line: Left: Wave- vector positions for different types of magnetic order in the extended BZ. The inner 
hexagon marks the first BZ. Middle: Examples for the fiow behavior of the leading susceptibility component. J2 — 0.1 shows 
a representative example of an RG fiow displaying order and J2 = 0.3 for a disordered phase. Right: Staggered and plaquet 
dimer response along the J3 = axis. The insets show the staggered and plaquet dimerization patterns. 



gion by increasing J2 , we observe the appearance of clean 
magnetic response peaks which we interpret as the onset 
of weak magnetic order. From the peak positions in k- 
space we divide this region into a collinear ordered phase 
for large J3 and a spiral ordered phase for smaller J3 
(Fig. [1]). Throughout parameter space, the spiral phase 
is very close to 120°-Neel order on both honeycomb sub- 
lattices except for a small region near the J2 axis where 
the wave vector deviates from commensurability (Fig.[2|. 
From the perspective of the degenerate classical spiral 
phase, this corresponds to a pinning of the order due to 
quantum fluctuations. A pronounced jump of the leading 
susceptibility wave- vector is seen as we cross the transi- 
tion between AFM and C-AFM order, pointing to a first 
order transition. The observations are consistent with 
the classical result except for the fact that the transi- 
tion between AFM and C-AFM order is shifted towards 
higher J2,c ~ 0.57 compared to jciassicai ^ q 5 
quantum corrections included in our calculation. 

We now focus on selected cuts through parameter 
space. We first investigate how the fiuct nation profile 
changes along the J3 = line (Fig. [2|. For small J2 
AFM order manifests itself in peaks at the corners of 
the extended BZ and a characteristic instability break- 
down of the fiow (lower line of Fig. [2| . As we increase 
J2, the AFM peaks rapidly decrease and from the dis- 
appearance of unstable fiow behavior we estimate the 
transition to be at J2 ~ 0.15. Inside the paramagnetic 
phase, such as at J2 = 0.3 depicted in Fig. |2j no clear 
peak structure is visible and the susceptibility fiow re- 
mains stable up to A ^ 0. Around J2 ~ 0.6 spiral order 
peaks emerge at wave vectors slightly shifted from the 
commensurate positions towards larger |k|. This feature 
is consistent with a large-S* expansion [M] which allows 



to select specific wave vectors out of a classical mani- 
fold of degenerate momenta. Upon further increasing J2 , 
the peak positions approach commensurability, i.e., the 
sublattices effectively decouple and individually exhibit 
120°-Neel order. During the fiow, these susceptibility 
peaks emerge at very small A-scales, giving us indication 
of a very weak magnetization. While the transition be- 
tween paramagnetism and spiral order is a bit smeared 
out at small J3, the onset of response peaks occurs more 
abruptly at larger J3. 

To resolve more information about the correlations in 
the disordered phase, we compute the staggered and pla- 
quet dimer response, i.e., the (dimensionless) factor of 
amplification of an external dimer-field perturbation ex- 
erted on the system. In doing so we can distinguish be- 
tween parameter regimes of different dimer fiuctuation 
strength. As we sweep through the paramagnetic phase 
at J3 = we find that the staggered dimer response 
is dominant for higher J2 while plaquet and staggered 
dimer response compete for lower J2 (Fig. [2|. The abso- 
lute dimer response amplitudes are smaller for lower J2. 
This is in qualitative agreement with VMC [16 as well as 
with ED studies [11]. It supports the view that if at all 
the system may form a spin liquid phase around this do- 
main which is also the parameter regime related to the 
honeycomb Hubbard model from a strong coupling ex- 
pansion. There, charge fiuctuations which are neglected 
in ^ may be sufficient to destroy the comparably low 
dimer ordering tendency. 

In addition, we investigate parameter lines varying J3 
for intermediate J2 through the disordered regime. The 
fiuctuation profiles and dimer responses for the param- 
agnetic regime at finite J3 are shown in Fig. |3] We see 
that at small J3 the fiuctuation profiles for J2 = 0.5 and 
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FIG. 3. (Color online) J3 sweep of the susceptibility for J2 = 0.5 (upper line) and J2 = 0.6 (lower line). As J3 is increased, 
the disordered phase enters the AFM phase for J2 = 0.5 and enters the C-AFM phase for J2 = 0.6. Right: Corresponding 
staggered dimer responses (plaquet dimer responses are much weaker and not shown). We find an increased propensity to 
staggered dimer order for larger J2. 



J2 = 0.6 are very similar but differ more with increas- 
ing J3 until eventually the J2 = 0.5 line leads into the 
AFM ordered phase while collinear order emerges on the 
J2 = 0.6 line. The paramagnetic phase shows rather 
complicated susceptibility profiles which have lost most 
signature of ordering fluctuations. A typical feature is 
the ringlike shape as seen e.g. at (72,^3) = (0.6,0.5). 
An intuitive reason for a quantum disordered phase is al- 
ready indicated from the classical limit where the point 
{J2iJ^) = (0.5, 0.5) is tricritical with three competing or- 
dering tendencies. From the dimer responses along the 
J2 = 0.5 and J2 = 0.6 line we find that rather indepen- 
dent of J3, staggered dimer ordering tendency is more 
efficiently established by larger J2. 

Note added. When this manuscript was completed we 
became aware of an independent work providing an anal- 
ysis of model ([T]) through joint mean field and exact diag- 
onalization techniques [28] . Several similar findings show 
a good correspondence of both approaches. 
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